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ABSTRACT 

We present a detailed survey of the dynamical structure of the phase space around the new 
moons of the Pluto-Charon system. The spatial elliptic restricted three-body problem was 
used as model and stability maps were created by chaos indicators. The orbital elements of 
the moons are in the stable domain both on the semimajor axis - eccentricity and - inclination 
spaces. The structures related to the 4:1 and 6:1 mean motion resonances are clearly visible 
on the maps. They do not contain the positions of the moons, confirming previous studies. We 
showed the possibility that Nix might be in the 4: 1 resonance if its argument of pericenter or 
longitude of node falls in a certain range. The results strongly suggest that Hydra is not in the 
6: 1 resonance for arbitrary values of the argument of pericenter or longitude of node. 
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1 INTRODUCTION 

In 1930 C. Tombaugh discovered Pluto, the ninth planet of the Solar 
system. From its discovery until 2006, Pluto was considered the 
Solar System's outmost planet. In the last two decades, however, 
many objects similar to Pluto were discovered in the outer solar 
system, notably the scattered disc object Eri s, which is 27% more 
massive than Pluto jSrown & SchaUeJl200 7^. On August 24, 2006 
the lAU defined the term "planet" for the first time. This definition 
excluded Pluto, which the lAU reclassified as a member of the new 
category of dwarf planets along with Eris and Ceres. 

Pluto's first moon, Charon was found bv lChristv & HarringtonI 

( Il978h . which greatly facilitated the study of Pluto, since that dis- 
covery made possible a more accurate determination of Pluto's 
mass. In Table [T] a chronology of mass ratios fi — mc/mp and 
mass parameters fi — mc/(mp + mc) is shown, where mp 
and mc are the masses of Pluto and Charon, respectively. Sev- 



Table 1. Values for the Charon/Pluto mass ratios and mass parameters. In 
the first six rows values deiived form the barycentric wobble are listed, 
while in the last two rows those calculated from orbital fits using the dis- 
covery of the small moons. 



Reference 


Mp 


X 10 




M = 4^ X 10 




Null et al. (1993) 


0.837 ± 


1.47e - 


2 


0.7723t;;235e- 


2 
2 


Young &Binzel (1994) 


1.566 ± 


3.50e - 


3 




3 
3 


Null & Owen (1996) 


1.240 ± 


S.OOe - 


3 


1 1032+^-2*'^'=- 


3 
3 


Tholen & Buie (1997) 


1.100 ± 


6.00e - 


2 




2 
2 


FoustetaL (1997) 


1.170 ± 


6.00e - 


3 


1.0474t^;^«3e- 


3 
3 


Olkin et al. (2003) 


1.220 ± 


S.OOe - 


3 


1.0873ie;3i0e- 


3 
3 


Buie et al. (2006) 


1.165 ± 


5.50e - 


3 




3 
3 


Tholen et al. (2008) 


1.166 ± 


6.90e - 


3 


l-0442lt;500e- 


3 
3 


Mean 


1.1183 






1.05537 





eral autnors have struggled to o btam tms quantit 
ments of the barycentric wobble ('NulletalJl993l: 



■Yo ung & Bi nzell 

1994; Null & Owen 1996; Tholen & Buiel l 19971 ; (Foust et al. 1991 
Olkin et alj|200 3). In addition the Pluto-Charon system is remark- 



able, since in the Solar system Charon is the largest moon relative 
to its primary, with the highest mass ratio of 0. 1 166 dTholen et al.l 
l2008h (hereafter referenced to as T08). A visual representation of 
the mass ratios and the corresponding mass parameters is displayed 
in Figure [T] Horizontal solid lines give the error bars of the mass 
ratios with a vertical tick in the center denoting the best-fit value de- 
termined by the authors. Below each solid line the corresponding 
mass parameter is plotted as dashed lines, where the vertical tick 
denotes the mass parameter computed from the best-fit value of the 
mass ratio. We note that this is not at the center, as it is evident 
from the third column of Table [T] The largest error bar is given for 
lTholen&Bui3 ( ll997l) since their observations were not optimized 
for the determination of the Charon/Pluto mass ratio. The mean of 



the mass ratios and ma ss paramete rs are given in the last row of 
Table [T] In the works of iBuie et al.l |2p06) (hereafter referenced to 
as B06) and T08 the masses were derived from two-body and four- 
body fits, respectively. From Table[T]one can see that the mass ratio 
has considerably evolved and reached by now a well established 
value. In the present work the dependence of the phase space on 
the mass parameter is studied. 

Subsequent searches for other moons around Pluto had been 
unsuccessful until mid May 2005, when two new moons were dis- 
covered ( Weaver et al. 2005, 2006). The discovery of Pluto's new 
moons. Nix (provisionally designated by S/2005 P2) and Hydra 
(provisionally designated by S/2005 PI) rendered the system even 
more interesting. They orbit the center of mass of the system, which 
is very close to the Pluto-Charon barycenter The preliminary or- 
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Tholen et ol. ( 200B) 

Buie e t qL ( 2006) 

OIkin et Ql. (2 003) 

Foust et Ql, ( 1997) 
_ J _ I 
Tholen & Buie (1997) 



Null & Owen ( 1996) 

Young & Binzel (1994) 



Null et Ql. (1995) 



■ Charon/Pluto moss ratio 
mass parameter 



Table 2. Physical values of the moons, where m is the mass, p is the den- 
sity, D is the diameter and a is the visual geometric albedo (T08). Square 
brackets indicate assumed quantities. 



0.06 0.08 0.10 0.12 0.14 0.16 0.18 

Charon/Pluto moss ratio ond mass porometer 

Figure 1. Best-fit mass ratio jl values for the Charon-Pluto system along 
with the coiTesponding mass pai'ameters fi. The horizontal bars span plus 
and ininus one standard deviation from the best-fit values. 



bits computed bv lWeaveret"aD(l2005h (hereafter referenced to as 
W05) were based on just two observations separated by only three 
days, considerably less than a full orbit of either moon. Because 
of these constraints unique orbits could not be calculated from the 
available data, but the measured positions were consistent with 
nearly circular orbits in the orbital plane of Charon. On this as- 
sumption, preliminary orbital solutions yield a = 64700 ± 850 km 
and P = 38.2 ± 0.8 days for Hydra, and a = 49400 ± 600 km and 
P = 25.5 ± 0.5 days for Nix. Two-body orbit solutions for Nix 
and Hydra were computed by B06 using images taken of the Pluto 
system during 2002-2003 with the Hubble Space Telescope. Their 
use of data derived from prediscovery observations that span sev- 
eral orbits of all the moons made possible to compute unrestricted 
fits to the orbits of Nix and Hydra. According to the results, the 
orbital periods of Nix and Hydra are close to the ratio of 4: 1 and 
6: 1 with that of Charon, respectively, indicating mean-motion reso- 
nance. In this paper two bodies are in mean-motion resonance when 
n/n' = p/ip + q), where n and n' are their mean motions, p and q 
are small integers, where q is the order of the resonance. The mean 
motions are measured in the sidereal system. 

Some of the major conclusions of B06 are as follows: (;) Nix 
and Hydra are in nearly circular orbits, with eccentricities of 0.0023 
and 0.0052, respectively, (ii) The orbits of the small moons are al- 
most coplanar with Charon's orbit. (Hi) The orbital periods of Nix 
and Hydra are nearly commensurate with the period of Charon, but 
differ significantly from the exact ratios of 4: 1 and 6: 1, respectively. 
They argue that perhaps there are no resonances acting between the 
bodies. We note that in B06 the eccentricity of Charon was assumed 
to be zero; they argue that the eccentricity of the orbit of Hydra is 
significantly nonzero, unlike the orbits of Charon and Nix, which 
are consistent with zero eccentricity. As we will demonstrate the 
eccentricity of Charon substantially influences the phase space of 
the Pluto-Charon system. 

Although the two-body orbit solutions are good enough to pro- 
vide satisfactory agreement with the 2002-2003 data, T08 studied 
whether the direct perturbations might be too strong to permit a suf- 
ficiently accurate extrapolation forward to the 2015 New Horizons 
spacecraft encounter with Pluto. Also, with adequate data, a formal 
four-body orbit solution should yield the mass for each member of 
the system. 

Their results of a four-body orbit solution for the Pluto system 
based on published observations have established Icr upper limits 
on the masses of Nix and Hydra. The masses and some other phys- 
ical properties of the moons are shown in Table |2l These limits 
place lower bounds on their geometric albedos and upper ones on 



Name 


m [kg] 


P 


D [km] 


Q(V) 


Pluto 


1.304 xl022 


2.06 


[2294] 


0.61 


Charon 


1.520 xlO^i 


1.63 


1212 


0.34 


Nix 


5.8 xlQi'' 


[1.63] 


88 


0.08 


Hydra 


3.2 xlO^'^ 


[1.63] 


72 


0.18 



their densities. The Charon/Pluto mass ratio was determined from 
the four-body orbit solution. The best-fit value of 0.1166 is almost 
identical to that of B06, since both team have used the same data. 

The work of T08 led to the conclusion that the orbits of 
Charon, Nix, and Hydra are not quite coplanar, and the latter 
moons' orbit planes precess around the system's invariable plane 
with periods of 5 and 15 years. The orbital eccentricities are 
nonzero but small for all three moons when measured in a barycen- 
tric reference frame. The results of T08 did not provide evidence on 
any mean motion resonances between the three moons confirming 
the study of B06. 

The eccentricity for Charon disagrees with the zero eccentric- 
ity published in B06. The erroneous value was the result of on the 
one hand using a software which was developed to fit highly ec- 
ce ntric orbit and on the other hand the use of the orbit published 
in iTholen & Bui3 jl997h . It should be noted that the magnitude of 
Charon's eccentricity is very similar to the values determined from 
Pluto surface models of Tholen & Buie ( 1997) (see Table III on 
page 251). In the case of the albedo model where only Pluto was 
considered, the computed eccentricity was 0.003 ± 0.0005 which 
agrees very well with the value of 0.0035 calculated in T08. 

To perform the four-body orbit solution all 22 parameters (18 
orbital elements and 4 masses) were fitted simultaneously in the 
work by T08. During the computations, the solution settled into 
many different local minima, therefore without computing every 
possible orbit, it can not be guaranteed to find the absolute min- 
imum. The currently available data are insufficient to produce a 
unique minimum. The possibility of stellar occultation opportuni- 
ties and new Hubble Space Telescope observations will likely im- 
prove and refine the orbital elements of these new moons. IVIuch 
more precise determination of the orbits and masses of Nix and 
Hydra will be possible as the New Horizons spacecraft approaches 
the Pluto system in 2015. 

Nagy, Siili & Erdi (2006) studied the phase space of the Pluto- 
Charon system in the framework of the spatial circular restricted 
problem. The moons were treated as test particles and their semi- 
major axes, eccentricities and inclinations were varied. The com- 
putations were repeated for different initial mean anomalies. Their 
results showed that the region inside « 42000 km is unstable, thus 
no moon could exist there. According to their results both moons 
reside in the stable region of the phase space of the Pluto-Charon 
system and the upper limit for the eccentricities of Hydra is 0.17, 
while it is 0.3 1 for Nix. In the semimajor axis - inclination plane the 
4:1 and 6:1 resonances are clearly visible above « 20° and « 35°, 
respectively. Unfortunately, these values could not serve to place a 
more stringent upper limit on these orbital elements. 

The m ain goal of t his pa per is to extend the previous inves- 
tigations of iNagy et al.l | |2006|) using the spatial elliptic restricted 
problem. In their work the eccentricity of Charon was zero and the 
mass parameter of the system was 0.130137. Since then both of 
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Table 3. Orbital parameters derived by W05, Keplerian fits by B06 and those from four-body orbits solution by T08 typesetted in bold. The 
values are valid at epoch JD 2452600.5, mean equator and equinox of J2000. The parameters for the orbit of Charon are relative to Pluto while 
the orbits of Nix and Hydra are relative to the center of mass of the Pluto-Charon system. The values in parentheses in the semimajor axis 
column are given in units of the semimajor axis of Charon (A=l). Square brackets indicate assumed quantities. 



Body 


a [km] [A] 


e 


i [deg] 


u [deg] 


n [deg] 


L [deg] 


M [deg] 


T [day] 


Charon (B06) 
Charon (T08) 


19571.4(1) 
19570.3 (1) 


0.0 
0.0035 


96.145 
96.168 


157.9 


223.046 
223.054 


257.946 
257.960 


257.946 
237.006 


6.3872304 
6.38720 


Nix (W05) 
Nix (B06) 
Nix (T08) 


49400 (2.524) 
48675 (2.487) 
49240 (2.516) 


[0.0] 
0.0023 
0.0119 


[96.145] 
96.18 
96.190 


352.86 
244.3 


223.14 
223.202 


7 

123.14 
122.7 


261. lA 
15.198 


25.5 
24.8562 
25.49 


Hydra (W05) 
Hydra (B06) 
Hydra (T08) 


64770 (3.309) 
64780 (3.310) 
65210 (3.332) 


[0.0] 
0.0052 
0.0078 


[96.145] 
96.36 
96.362 


336.927 
45.4 


223.173 
223.077 


7 

322.71 
322.4 


122.61 
53.923 


38.2 
38.2065 
38.85 



these values have been updated and also the orbital elements of 
Nix and Hydra have been refined as it can be seen from Table [3] 
The relevant changes are the following: (;) the semimajor axes of 
both moons have been changed, (li) Nix's eccentricity is signifi- 
cantly different from zero, and Hydra's eccentricity has increased 
by a factor of 150% (Hi) the argument of periapses have completely 
different values. The most important change is the transition from 
the circular to the elliptic case. In the circular restricted problem 
there exists a first integral of motion, i.e. the Jacobi-integral, which 
does not exist when introducing the eccentricity. Due to the disap- 
pearance of the Jacobi integral the phase space structure changes 
qualitatively. 

With that said, the orbital elements tabulated in Table [3] are 
subject to change as new observations become available. Using sta- 
bility maps computed in advance for a large set of orbital param- 
eters has the advantage that the stability properties of the moons 
can be easily established when the orbital parameters are modified. 
Moreover the error bars can be readily combined with the stability 
maps giving extra information about the motion. 

In Section 2 we describe the investigated model and give the 
initial conditions used in the integrations. The applied numerical 
methods are briefly explained in Section 3. The results are pre- 
sented in Section 4. Section 5 is devoted to the conclusions. 



2 IVIODEL AND INITIAL CONDITIONS 

Since the orbital radii of the moons are much smaller than the Hill 
radius (~ 8.0 x 10* km) of the Pluto-Charon system, the moons 
are deep in Pluto's gr avitational well, the perturbations by the Sun 
can be ignored, as did lLee & Pe ale (2006) and T08. 

To study the structure of the phase space of the Pluto-Charon 
system we applied the model of the spatial elliptic restricted three- 
body problem. We integrated the dimensionless equations of mo- 
tion. An obvious advantage of using such equations is that the re- 
sults are independent of the exact value of the semimajor axis of 
Charon. The unit of length was chosen such that the separation of 
Pluto and Charon (the primaries) is unity, i.e. the semimajor axis of 
Charon oi = 1 A in all computations. Let the unit of mass be the 
sum of the primaries, i.e. rap + mc = 1, where mp and mc are 
the masses of Pluto and Charon, respectively. Let the unit of time 
be chosen such that fc'^(mp+mc) ~ 1, where fc is the gravitational 
constant. The orbital plane of the primaries was used as reference 
plane, in which the line connecting the primaries at t = defines 
a reference i-axis. The longitude of the ascending node is the 



angle between the line of nodes and the x axis. The argument of 
pericenter ua is measured between the line of nodes and the radius 
vector of the pericenter. The values of f2 = 0, = and the mean 
anomaly M = Q places the test particle on the x axis at a distance 
of a(l — e) from the barycenter, where a is the semimajor axis and 
e is the eccentricity of the test particle. 

The initial orbital elements of Charon are given in Table [3] 
where those of T08 were used. Hereafter the orbital elements of 
Charon will be denoted by subscript 1, i.e. ei is the notation for 
Charon's eccentricity. 

Though Nix and Hydra are almost coplanar with Charon, still 
we study the problem more generally by considering the effect of 
non-zero inclinations on the orbital stability. The mass parameter 
^ = 0.104424 was chosen according to the value of 0.1166 pub- 
lished in T08 (see Table[T). 

To examine the phase space in the vicinity of Nix and Hydra 
separately we varied the initial orbital elements of the test particles. 
Stability maps were created for the (a — e) and (a — i) orbital ele- 
ment space for both moons for different eccentricity e, inclination 
i, argument of pericenter uj and longitude of node Q.. The values 
of the mean anomalies M were kept constant in all the simulations 
and are given in Table [3] To compute the (a — e) and (a — i) stabil- 
ity maps the semimajor axis was varied with a stepsize of lO^"' in 
the intervals which are given in Table|4] For each computed (a — i) 
map the inclination was changed from to 45° with Ai = 1° and 
for each (a — e) map the eccentricity varied from to 0.4 with 
Ae = 10~^ 

The (a — e) stability maps were computed for different pa- 
rameters of the moons. For each {[i]k, and {[i\k, pair an 
(a — e) stability map was computed, where 

[i]^ = kM', fc = 0, ...,4, 

= 22.°9 + /Atj', [J7],=ZAn', Z = 0, ...,7, 

where Ai' = 10° , Atj' = Afi' = 45°. Here we used the [] notion 
to distinct the variables from those of Charon. Computation of each 
combination of the (fc, I) integer pair results in 2 x 40 stability maps 
per moons. 

Similarly the [a — i) stability maps were computed for each 

([e]fc, [u)]i) and ([e]fe, pair, where 

[e]^ = kl\e, fc = 0, ...,4. (1) 

where Ae' = 0.1. The details are summarized in Table|4] 

In total more than 14 million orbits were calculated, which re- 
sulted in 2 X 80 = 160 stability maps for each moon. Due to the 
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Table 4. The initial orbital elements for the test particles. In the upper part 
the intervals of the semimajor axes along with the stepsizes A are listed for 
the stability maps (a — e) and (a — i). The lower part shows the intervals 
(I) for the orbital elements [e] ^ , [oj] ; and [Q] i and the respective stepsizes 
(A'). 



Map 


a [A] 


e 


i [deg] 




Nix 


[2.40,2.64] 


[0,0.4] 


[0,45] 




Hydra 


[3.22,3.38] 


[0,0.4] 


[0,45] 




A 


10-^ 


10-3 


1 






[e]fe 


[i]k [deg] 


[deg] 


Ml [deg] 


I 


[0,0.4] 


[0,40] 


[22.9;337.9] 


[0;315] 


A' 


0.1 


10 


45 


45 



very small stepsize in a and e, each (a — e) stability map corre- 
sponds to more than 6 x 10'' initial conditions, thus providing a 
very fine resolution. 

The above orbital elements refer to a barycentric reference 
frame, where the mass of the barycenter is mp + mc- By the usual 
procedure we calculated the barycentric coordinates and velocities 
of the test particle and then transformed them to a reference frame 
with Pluto in the origin. In the numerical integrations we used the 
latter coordinates and velocities. 

For the integration of the system we applied an efficient 
variable-timestep algorithm, known as Bulirsch-Stoer integration 
method. The most important feature of this algorithm for simula- 
tions is that it is capable of keeping an upper bound on the local 
errors introduced due to taking finite timesteps by adaptively reduc- 
ing the step size when interactions between the particles increase in 
strength. The parameter e which controls the accuracy of the inte- 
gration was set to e = W~^°. The orbits were integrated for 10'' 
Charon's period (hereafter Tc). 



3 IVIETHODS 

To compute the stability maps, the method of the maximum eccen- 
tricity (ME), the Lyapunov characteristic indicator (LCI) and the 
relative Lyapunov indicator (RLI) were used as tools for stability 
investigations of the massless bodies representing the small moons. 

The IVIE method uses as an indication of stability a straight- 
forward check based on the eccentricity. This action-like variable 
shows the probability of orbital crossing and close encounter of 
two bodies and therefore its value provides information on the sta- 
bility of orbits. This simple check was already used in several sta- 
bility investigations, and was found to be a powerful indicator of 
the stability character of orbits l lDvorak et alj2003l : ISuli et alj2005l : 
iNagv et alj2006l) . In this work we define IVIE as follows 

ME = max (e). 

t£lO,10-i]Tc 

As a complementary tool, we computed also the LCI, a well- 
known chaos indicator. The LCI is the finite time approximation of 
the largest Lyap unov exponent, which is described in detail in e.g. 
iFroeschlg il984) . The definition of the LCI is given by 

LCI(i,xo,eo) = ilog||C(t)||, 

where xo is the initial condition of the orbit and ^(t) is the solution 
of the first order variational equations. The function LCl{t, a;o, ^o) 
measures the mean rate of divergence of the orbits. 



In lSandor et all ( |2000|.|2004|) the difference between the LCIs 
of two neighbouring orbits was introduced: 

AL(xo, Co, A) = |LCI(xo, Co, t) - LCI(a;o + A, Co, t)\, 

where A is the distance between two close orbits. This quantity 
measures the fluctuations of the curve of the LCI. To smooth the 
time evolution of AL its time average is computed, which is the 
definition of the RLI: 

RLI(xo,eo,A) = i^|AL(:co,Co,A)|, 

The definition contains A as free parameter, which must be chosen 
small enough to refle ct the local properties of th e flow in the phase 
space. It was shown, dSandor et al.l l l2000ll2004h ) that the choice of 
this parameter in a quite large interval (A G [10"''', 10"*^]) does 
not modify essentially the behaviour of the RLI. This method is ex- 
tremely fast to determine the ordered or chaotic nature of orbits and 
the method is sensitive for sho wing the resonant structure on sta- 
bility maps dSandor et ai]|2006h . In our computations A = 10 
was used. 

The three methods are not equivalent, however they complete 
each other. For example, the ME of the Earth is small, indicat- 
ing stability, although we know from numerical experiments that 
in fact the Earth is moving on a chaotic trajectory with a small 
but nonzero Lyapunov exponent. Therefore the ME detects macro- 
scopic instability (which may even result in an escape from the 
system), whereas the RLI and the LCI are capable to indicate mi- 
croscopic instability. 

According to notation. (??) the stability map is prepared in 
such a way, that for each set of the described initial conditions the 
LCI, RLI and ME value of the corresponding orbit of the test par- 
ticle is computed and plotted in the (a — e) or (a — i) parameter 
plane. 



4 RESULTS 

In this Section the ME_ stability maps are presented (Figs |2}j7]l, 
where ME_ is defined as follows: 

f 1 ifME = l 

ME_ = < 

ME — eo otherwise 

where eo is the initial eccentricity of the test particle. 'We used the 
above definition in order to display the maximum change in e in 
the course of the integration and to be able to compare orbits with 
different eo. For each (a, e) or (a, i) point the ME_ is computed 
and depending on its value a greyscale dot is plotted. Darker shades 
correspond to very low values of the ME_ and regular behaviour 
of the test particle, while lighter shades indicate larger ME_. In 
order to better visualize the structures in the stability maps the scale 
for the (a — e) figures extends from zero to 0.2 for Nix and 0.1 
for Hydra, for the (a — i) maps to 0.2 and 0.05, respectively. In 
the (a — e) maps the thick black solid curve in the white upper 
region is a contour line where ME is approximately equals to 0.9 
and marks the upper boundary of stable region. Above the curve 
ME is practically 1.0 while below it ME rapidly decreases. The 
thin dashed contour line right below the thick one corresponds to 
ME=0.4 and in all regions below this curve the ME_ is less than 
0.2 and 0.1 for Nix and Hydra, respectively. In the (a — i) maps 
ME is less than 0.4 in all cases. 

On the figures white solid curves are contour lines: along these 
curves the ME_ has a constant value. These values were chosen 
in such a way that the curves draw the approximate boundary of 
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the most prominent structure and to make comparison possible. 
The white numbers at the upper end of these curves are the cor- 
responding uj or Q, values. iMurrav & Dermot 3 (Il999h derived the 
maximum libration zone of resonances as a function of the semi- 
major axis and eccentricity using an analytical model based on the 
circular restricted three-body problem. These zones for g ^ 2 reso- 
nances have a V-shape on the (a — e) plane. The V-shapes are well 
represented by the white contour curves. Using a different, but rea- 
sonable ME- value would provide a similar V-shape. The chosen 
values to plot the contour curves are based on several tests. 

The stability maps based on the LCI and RLI values are essen- 
tially the same as those based on ME_, therefore these maps are 
not presented. 



4.1 The mass parameter 

The mass parameter used by iNagv et al.l j2006h differs by 25% 
from the present best-fit value therefore we studied the system with 
different ^ values. For the sake of comparison four runs were per- 
formed, one for /i = 0.130137 and one for ^ = 0.104424 in the 
(a — e) phase space in the vicinity of both Hydra and Nix. The 
stability map in the vicinity of Nix for /i = 0.104424 and ei = 
is shown in the upper panel of Figure |2] The phase space struc- 
tures for the two ^ parameters (the one for ^ = 0.130137 is not 
shown) are almost identical, they are only shifted along the hori- 
zontal axis: structures for the lower mass parameter are closer to 
Pluto. This is a consequence of the changes in the coordinates of 
the barycenter. This shift amounts approximately to 0.005 A. The 
ME=0.9 contour line is almost the same although in the boundary 
zone between the unstable and stable region the maximum differ- 
ence in ME- computed using the two mass parameters can reach 
values as high as ~ 1. This large difference is a natural conse- 
quence of the chaotic nature of orbits emanating from this region. 
Using longer integration time would cease these large differences. 
The conclusion is that in this mass parameter region of ~ 0.1, 25% 
change in shifts the phase space structures along the horizontal 
axis. This shift could be important in the close vicinity of mean 
motion resonances, like in the case of Hydra where the lower end 
of the 6:1 resonance for /i = 0.130137 is very close to T08 but it 
shifts away for ^ = 0.104424 (see Fig.O. 



4.2 Charon's eccentricity 

Since the Jacobi-integral does not exist in the elliptic problem 
therefore the phase space structure is qualitatively distinct in the 
two models. In order to visualize the effect of the disappearance of 
the Jacobi-integral, the circular and the elliptic cases for Nix are 
plotted in Fig. |2] and for Hydra in Fig. [3] where the circular case is 
shown in the upper the elliptic case in the lower panel. For Nix two 
major differences comparing with the circular case are the follow- 
ing: (0 the unstable zone is much larger in the elliptic case begin- 
ning from e ~ 0.17 onward and (n) the center of the 4:1 resonance 
is shifted from a ^ 2.544 A to a ^ 2.564 A. The shape of the 
resonance did not change significantly but the fine inner texture is 
distinct. For Hydra the upper unstable zone extends too, the cen- 
ter of the resonance moves closer to Pluto and its shape changes 
completely. In both cases the resonance becomes stronger for low 
eccentricities in the elliptic case. 
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Figure 2. The (a — e) stability map for Nix with /i = 0.104424. In the 
upper panel the circular (ei = 0), in the lower panel the elliptic case 
(ei = 0.0035) is displayed. Orbits emanating above the upper solid thick 
black curve are unstable orbits with ME^ 0.9. The plus signs con'espond 
to the orbital elements published by B06 and T08. Along the contour curves 
ME_ = 0.04. 



4.3 The (a — e) stability maps 

The (a — e) stability maps are shown in the lower panels of Fig.|2] 
for Nix and in Fig.|3]for Hydra. The orbital elements of the moons 
were taken from Table[3]row T08. Each of the two figures is dom- 
inated by a V-shaped gray structure, corresponding to the 4: 1 and 
6: 1 mean motion resonances between Charon and Nix and Hydra, 
respectively. These resonances can represent either ordered (stable 
for infinite time), or weakly chaotic (which may become unstable 
after very long time) behaviour. The location of the plus signs cor- 
responds to the orbital elements published by B06 and T08. The 
size of the plus sign is not proportional to the error bars of the a 
and e values. The white curve is a contour line corresponding to 
ME_ value of 0.04 for Nix and 0.025 for Hydra. 

From Fig.|2]it is clearly visible that the present position of Nix 
is not in the 4: 1 mean motion resonance, although it is closer than 
the position of B06. Our results are in agreement with B06 and T08 
since none of these works have identified any resonant arguments 
between Nix and Charon. We note that there is a sharp jump in 
the ME- values along the left separatrix while the ME_ decreases 
smoothly as a function of a as the right separatrix is approached. 

Fig. [3] shows the stability map for Hydra. The V-shape 6: 1 
resonance is clearly visible and its lower peak is very close to 
the present position of Hydra. However, the ME_ values are very 
small, the white contour curve indicates the region where ME_ is 
above 0.025. The ME- in the close vicinity of Hydra is less than 
0.02 indicating that this high order resonance is very weak in the 
close proximity of Hydra. Again this result confirms those of B06 
and T08. 

In the next step we investigated the orbits systematically by 
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Figure 3. The (a - e) stability map for Hydra with ^ = 0.104424. The V- 
shape structure corresponds to the 6:1 resonance. Along the contour curves 
ME_ = 0.025. See explanation of Fig.|2] 

changing the initial orbital elements of the test particle as described 
in Section 3 (see also TableQ. In Fig.|4]the results are summarized 
for the 8 values of to. These "maximum" stability maps were ob- 
tained as follows: for each (a, e) point we plotted the maximum 
selected from the ME_ values computed for the 8 different u. The 
white curves on the figures approximately denote the location of 
the resonances as they appear on individual stability maps, like in 
Figs|2]and[3] 

In Fig. |4] the upper panel shows the parameter space around 
Nix where the contour lines belong to ME_ = 0.04. In the case 
of u — 157.° 9 Nix is inside the 4:1 resonance, and between the 
curves the ME_ values are significantly higher than for the other 
values of oj. In this case the pericentrum of Nix is the same as that of 
Charon. The detailed study of this configuration will be carried out 
in a later work. The former B06 position is obviously outside of the 
resonance confirming the findings of B06. The location, size and 
shape of the resonance strongly depends on the u of Nix. From Fig. 
|4]it is evident that the effect of the 4:1 resonance is not negligible 
for even very small eccentricities. 

In Fig.|4]the lower panel displays the parameter space around 
Hydra where the contour lines belong to ME_ — 0.025. The 
ME=0.9 curve depends approximately linearly on a, starting from 
e = 0.3 at a = 3.26 A and ending at e = 0.33 at a = 3.38 A. 
The 6:1 mean motion resonance is visible for all ui, but it can not 
be traced down below e « 0.05 for any value of the pericenter. The 
resonance is the strongest for u — 67° in the sense that its effect 
is visible for the lowest e ~ 0.05. The present position of Hydra 
is closest to the resonance when cj — 22°: if the upper limit of 
the scale is reduced to 0.03 than the resonance lowest peak is just 
above the T08 plus sign. As it is obvious from Fig.|4]the 6: 1 mean 
motion resonance is very weak in the vicinity of the moon, no sign 
of it can be observed around T08 or B06. 

The Lu published in T08 differs by 108° from that of B06 in 



the case of Nix, and by 68° in the case of Hydra (see Table[3j. Since 
the eccentricity is small for all three moons the value of lj is hard to 
establish accurately. New and more accurate orbital solution based 
on data yielded by future observations of the system could result in 
such Lj that would place Nix in the resonance. In the case of Hydra 
this is questionable since the 6:1 resonance is very weak for low 
eccentricities. The same applies to the nodes because the moons 
are nearly coplanar, however the value of Q changed only a few 
arcminutes. 

We note that the switch from 2D to 3D in the case of small 
inclination (i ^ 0.°1) does practically not induce any observ- 
able changes in the (a — e) parameter space. Computations were 
performed for zero inclination and the comparison of the stability 
maps have revealed only tiny differences. 

We briefly give a summary for higher inclination in the follow- 
ing. Nix's maximum stability map for i = 10° indicates stronger 
effects of the 4: 1 resonance and the expansion of the unstable zone 
which starts already from e = 0.1. The resonances are wider and 
closer to Pluto and to the T08 position, which is now inside more 
resonance structures indicating the potential of active resonance. 
Increasing the i up to 20° further shrinks the stable domain which 
ends at e ~ 0.06. The resonance moves even closer to Pluto and 
gain in power. The maximum stability maps for i — 30° and 40° 
show that the lower border of the unstable zone stays around 0.06 
and the resonance structures move closer to the barycenter. 

Hydra's maximum stability map for i = 10° is similar to Fig. 
|4l the ME=0.9 contour curve is practically the same as well as the 
resonant structures which have become slightly more prominent. 
For i ^ 20° the resonant structures can be followed down to the a- 
axis and like in the case of Nix they gain in power and shift closer to 
Pluto: in the i — 40° case the T08 position is outside the resonance. 
The upper limit of the stable zone is e ~ 0.28 at i = 40°: large 
increase in the inclination did not significantly reduce the size of 
the stable domain around the moon. 

For both moons the "maximum" stability maps for the 8 values 
of Q were obtained on the analogy of Fig.|4] For i = jn and i — 
iH the results agree quite well with the ones obtained by ui: each 
structure on Fig. |4] has a counterpart with JIn ~ (^n — 22°, and 
SIh ~ <jJh — 180°. The explanation of these simple relationships 
is the following: in the limit as i — > the orbital plane coincides 
with the reference plane and we have w = lj + i}, where -oo is the 
longitude of the pericenter. In the case of Nix from the equation 

zu^ = vjfn, 

where zuu, = + SIn and zuq = cjn + it follows that 
Qi — uii — 21°, where i = 1, . . . , 8. The equation yields for Hydra 
Qi^uj^- 182°. 

As a consequence the stability map belonging to SI = 135° 
indicates that Nix might be in the 4:1 resonance. In the case of 
nearly coplanar orbits even small details are similar on both figures. 
The above relationships are also true for i — 10° for the dominant 
features but the details are already distinct. The results disagree for 
higher inclination. 

4.4 The (a — i) stability maps 

The (a — i) stability maps are shown in Fig. |5]for Nix and in Fig. 
[6] for Hydra. The orbital elements of the moons were taken from 
Table |3] row T08. In the case of Nix the contour line belongs to 
ME- = 0.08; we note that this value is the double than those 
belonging to the (a — e) maps. In accordance with Fig.|2]the results 
displayed in Fig.|5]confirms that the present position of Nix is not 
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Figure 4. The (a — e) stability map for 8 uj values for Nix (upper) and Hydra (lower). The white curves correspond to contoiu" lines indicating the shape 
of the resonances as they would appear on individual maps for the ui displayed at the upper end of the curves. In the upper panel along the contour curves 
ME_ = 0.04. In the case of Hydra along the contour curves ME_ = 0.025. 



in the 4:1 mean motion resonance, but in a region where ME_ 
is very small indicating stable motion. The remarkable feature of 
the resonance is a vertical, approximately 0.008 A wide and 10 
degree high column at o = 2.564 A. Let us note that the shape of 
the resonance is reminiscent of a bee's abdomen with a sting. As 



before there is a sharp jump in the ME_ values along the left side 
and ME- decreases smoothly towards the right separatrix. 

For Hydra the top of the scale of ME_ is 0.05 and the white 
contour line belongs to ME_ = 0.017. The 6:1 mean motion res- 
onance structure is observable only for i ^ 20° where the maxi- 
mum of ME- is less than 0.05 indicating that the resonance is very 
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Figure 5. The (a — i) stability map for Nix witli initial conditions taken 
from T08. Along the contour curves ME_ = 0.08. 
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Figure 6. The (a — i) stability map for Hydra with initial conditions taken 
from T08. Along the contour curves ME_ = 0.017. 



weak. In accordance with Fig.[3]Hydra is not in the 6:1 mean mo- 
tion resonance. The two regions enclosed by the white curves on 
the right and left side of the resonance contain values that are less 
than 0.017. It is interesting that two local minima appear at i ~ 40° 
on both sides of the resonant structure. 

In Fig. |7] the maximum stability maps for Q for both moons 
are presented: the upper panel for Nix, the lower one for Hydra. 
The Q = 135° map shows that T08 of Nix is in the "sting" of 
the 4:1 resonance while all others are at larger semimajor axes 
(Fig. [7] upper panel). In the case of Hydra there is no sign of 
any active resonance in the vicinity of T08, all the observable ef- 
fects are above 10 degrees. The most prominent features belong to 
n = 180° 225° 45° and to 0°. 

We have performed runs for e = in order to estimate the 
effect of the small eccentricities. It turned out that no observable 
difference could be seen when comparing it with the eccentric case. 

The maximum stability map for Nix with initial e = 0.1 pre- 
dicts a very narrow stable zone between and i ^ 2° at the posi- 
tion of T08. According to the upper panel of Fig.|4]the motion with 
initial e > 0.15 is completely unstable. 

The 6:1 resonance with initial e = 0.1 is stronger and the T08 
position of Hydra is very close to it. The investigated domain is sta- 
ble up to 45°, the ME- values rise above 0.1 beyond 40° but no 
escape occurs. Interestingly the parameter space is more stable for 
higher eccentricity (e = 0.2) and unstable motion occurs only for 
e ^ 0.3, although in accordance with the lower panel of Fig.|4]or- 
bits with a > 3.28 A are stable if the inclination is small {i ^ 10° ) . 
In the case of e = 0.4 shear instability characterize the system. 



5 CONCLUSIONS 

The phase space around the two small moons. Nix and Hydra of 
the Pluto-Charon system was studied in detail using the framework 
of the elliptic restricted three-body problem. As initial conditions 
the orbital elements of T08 were used and were varied to map the 
moons' phase space. We studied the effect of increasing the mass 
parameter and found out that 25% change only marginally influ- 
ences the results: the most significant effect is a shift of the struc- 
tures along the horizontal axis. This might be important in the cases 
of resonances. 

The former work of iNagy et all ilOOdi) used the circular prob- 
lem, but present observations favor elliptic orbit of Charon (T08). 
To compare the elliptic and the circular case computations were 
performed in both models. Several important differences were ob- 
served: (0 the unstable zone is much larger in the elliptic case, (ii) 
the center of the 4:1 resonance is shifted from a « 2.544 A to 
a ~ 2.564 A, although its shape did not change significantly, (Hi) 
the 6: 1 resonance is shifted too and it has a completely different 
shape and (iv) the resonances become stronger for low eccentrici- 
ties in the elliptic case. 

The present positions of the moons (denoted by T08 on the 
figures) are in stable regions both on the (a — e) and (a — i) param- 
eter spaces. On the stability maps the structures related to the 4: 1 
and the 6: 1 mean motion resonances are clearly visible but none of 
them contains any of the moons. These results are in line with those 
of B06andT08. 

"Maximum" (a — e) stability maps were created for both 
moons. The map for Nix shows the possibility of active resonance 
in the case ofu; — 157.° 9 or f2 = 135° . For Hydra the "maximum" 
stability map did not show this possibility, since the 6: 1 resonance 
is very weak in the vicinity of the moon. It was also showd that the 
planar and spatial cases are almost identical, when the inclinations 
are very small, like those of the moons. 

Analogous maps were presented for the (a — i) plane which 
confirmed our previous findings. Interestingly in the case of Nix the 
contour line belongs to 0.08, which is the double of that belonging 
to the (a — e) maps. 
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